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Abstract 

A methodology for accurate and efficient simula- 
tion of unsteady, compressible flows is presented. The 
cornerstones of the methodology are a special discretiza- 
tion of the Navier-Stokes equations on structured body- 
fitted grid systems and an efficient solution-adaptive 
mesh refinement technique for structured grids. The 
discretization employs an explicit multidimensional up- 
wind scheme for the inviscid fluxes and an implicit 
treatment of the viscous terms. The mesh refinement 
technique is based on the AMR algorithm of Berger and 
Colella ( J. Comp. Phys., Vol. 82, pp. 64~84> 1989). In 
this approach, cells on each level of refinement are or- 
ganized into a small number of topologically rectangu- 
lar blocks, each containing several thousand cells. The 
small number of blocks leads to small overhead in man- 
aging data, while their size and regular topology means 
that a high degree of optimization can be achieved on 
computers with vector processors. 


Introduction 

Many flows of interest to scientists and engineers 
are fundamentally unsteady in character. Yet, when 
computational methods are used to analyze such flows, 
the effects of unsteadiness are often neglected and the 
flows are simply modeled as steady. The use of Com- 
putational Fluid Dynamics (CFD) in the design of gas 
turbine engines is a prime example of this practice. The 
flow of gases through multi-stage compressors and tur- 
bines is typically modeled, stage by stage, as steady. 
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Any unsteadiness due to passage of wakes from up- 
stream blade rows is ignored. At best, the effects of 
unsteadiness is accounted for through empirical cor- 
relations (see e.g., Ref. 1). The reason behind this 
“pragmatic” approach is that with current computer 
resources and CFD algorithms, accurate simulations of 
unsteady viscous flows in turbomachinery are forbid- 
dingly expensive. CFD algorithms that minimize the 
cost of such simulations while guaranteeing sufficient 
accuracy are needed. 

For simulating steady and unsteady inviscid flows, 
techniques that combine high resolution shock captur- 
ing schemes with unstructured grid systems and solu- 
tion adaptive mesh refinement have proven quite suc- 
cessful in recent years (see, e.g., Ref. 2-5). Most of 
these techniques tessellate the flow domain using either 
fully unstructured grid systems (triangles and tetrahe- 
dra) or a patchwork of non-boundary-conforming Carte- 
sian grids. The attractive feature of such grid systems 
is that they can be generated quickly and autonomously 
for geometries of arbitrary complexity, significantly re- 
ducing the man-hours needed for problem setup. 

For viscous flows, these same techniques have not 
yet brought the same success. Several factors con- 
tribute to this lack of success. First, adequately resolv- 
ing boundary layers in high Reynolds number flows has 
proven difficult. 6 ’ 7 Second, the non-trivial data struc- 
tures required for unstructured grids do not lend them- 
selves to high level of optimization on computers with 
vector processors. This can be a particularly serious 
drawback in simulations of viscous flows which require 
far greater resolution of the flow field than do inviscid 
flows and, hence, greater computational efficiency. The 
difficulty is often exacerbated by a significantly greater 
memory usage of computer codes employing unstruc- 
tured grids compared to codes employing structured, 
body-fitted grid systems. Finally, performance of im- 
plicit discretizations on unstructured grids has not been 
satisfactory. 6 While this is not a serious limitation for 
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inviscid flows, it can be a major limitation in simu- 
lations of viscous flows due to cell-Reynolds number 
constraint on the allowable time step size of explicit 
schemes. 

Compared to unstructured and non-boundary-con- 
forming Cartesian grids, structured body-fitted grid sys- 
tems offer many advantages related to accuracy and 
efficiency in simulations of viscous flows. Specifically, 
boundary layers can be resolved very efficiently with 
structured body-fitted grids. Also, the regularity of the 
grid system admits simple array data structures, facil- 
itating optimization and vectorization of the computer 
codes. The regularity of the grids likewise facilitates 
application of implicit discretizations. Finally, genera- 
tion of structured grids is becoming easier than ever be- 
fore. Recent developments in techniques for generating 
high quality, structured, body-fitted, multiblock grid 
systems ®" 12 indicate that same level of autonomy can 
be achieved in generation of structured grids for compli- 
cated geometries as in generation of unstructured grids. 

To minimize cost of simulations that use struc- 
tured grid systems and to enhance acuracy of the com- 
puted results, techniques need to be developed for per- 
forming local solution-adaptive refinement in structured 
grids. Such techniques must avoid the complicated data 
structures on the cell level that are used in unstructured 
grids, or else many of the important advantages articu- 
lated above will be lost. One grid refinement technique, 
generally applied to Cartesian grids, appears exception- 
ally suitable for structured grids. This technique is the 
AMR (Adaptive Mesh Refinement) algorithm of Berger 
and Colella . 13 The AMR algorithm takes advantage of 
the fact that, on a given grid system, cells that require 
refinement come in clusters — a whole area needs refine- 
ment rather than single cells. The AMR algorithm or- 
ganizes these clusters of cells into a small number of 
topologically rectangular blocks, each containing a few 
thousand cells. Thus, simple and efficient array data 
structures can be used for storage of data in each grid 
block and a block-structured grid can be maintained on 
every level of refinement. The relatively small number 
of blocks leads to small overhead in manipulation of 
data during computations. 

The objective of this research is to develop an 
efficient methodology for performing accurate simula- 
tions of unsteady viscous flows that capitalizes on the 
advantages of structured grid systems. Cornerstones 
of the new methodology include a special discretiza- 
tion of the compressible Navier-Stokes equations that 
is designed especially for accuracy and efficiency, and a 
solution-adaptive mesh refinement technique for struc- 
tured, curvilinear grid systems that is based on the 


AMR algorithm of Berger and Colella . 13 The discretiza- 
tion employed is second order accurate in both space 
and time. It combines the explicit multi-dimensional 
upwind scheme of Colella 14 for the inviscid fluxes with 
an implicit scheme for the diffusion terms. The multi- 
dimensional upwind scheme has proven highly accu- 
rate in computations of inviscid flows , 14 ’ 5 while the 
implicit treatment of the diffusion terms significantly 
improves the robustness and computational efficiency 
of the scheme by eliminating the cell Reynolds num- 
ber constraint on the time step size that results from 
explicit treatment of the viscous terms. The use of solu- 
tion adaptive mesh refinement enhances both accuracy 
and efficiency of the scheme by providing high mesh res- 
olution in regions where it is required, but only where 
it is required. 

This paper focuses on only a few main aspects 
of this ongoing research. These aspects are the dis- 
cretization of the governing equations, the extension 
of the AMR algorithm to structured, body fitted grid 
systems, and implementation of the AMR algorithm 
using a mixed language programming (C++ and FOR- 
TRAN). Some sample results are shown at the end of 
the paper. 

Governing Equations and Discretization 

The governing equations employed are the com- 
pressible Navier-Stokes equations (see e.g Ref. 15). 
The fluid is taken to be a calorically and thermally 
perfect gas. Transformed to general curvilinear coor- 
dinates (£, 17 ) and cast in conservation law form, the 
governing equations can be written as 

d(JU) dF< dF« dFS dF? _ n 

dt + d£ + di) d£ di j 1 j 

where U = (p pu pv pe) T is the vector of conserved 
variables: density, momentum per unit volume in the 
x— and y — coordinate directions, and total energy per 
unit volume, respectively, and J is the transformation 
Jacobian. F ^ and F v represent the inviscid fluxes in 
the £ and 77 directions, respectively, whereas F$ and F% 
represent the viscous fluxes. Temperature is related to 
the remaining variables through 

r=( 7 -l)(e-i(« 2 + * 2 )) (2) 

For details of the fluxes see e.g. } Ref. 16. 

The governing equations are discretized using a 
hybrid explicit-implicit, cell-centered, finite volume dis- 
cretization on curvilinear, structured grids. The dis- 
cretized equations for cell (i,j) are written as follows: 
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The operator L can be written as 


up+ l - un 
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+ LW n ) 


-\{LUU n ) + L tt ij (U n+1 )) =0 


(3) 


where is the area of cell (t,;), L c {j represents the 
discretization of the convection terms in Eq. (1) and 
£V. represents the discretization of the viscous terms. 
The operator L is derived using a straightforward 
extension of the multi-dimensional upwind method of 
Colella 14 to viscous flows. In this method, 2% is written 
as follows: 




where Ao< + i/ 2 ,i is the area °f cell f ace (* + 1/2, j) (be- 
tween cells (i,j) and (i + 1,;)), Aa^j+ 1 / 2 , is the area of 
cell face ( i,j + 1/2) (between cells (i,j) and (t, ;' + 1)), 
and and E^Y /2 are the numerical fluxes nor- 

mal to those faces. To illustrate how the fluxes are eval- 
uated, consider cell face (1 + 1/2, j). The flux 
is evaluated using an approximate Riemann solver and 
can be written as: 


F ?*',Z = fiO O “<«/«) <*> 

where subscripts “Z” and “72” indicate “left” and “right” 
state for the Riemann solver, n-i±i/ 2 t j ls th e unit normal 
to the cell face at (i ± 1/2, j) and 


L v j{U ) = [ L$(V) ) (7) 

V Lij (T) + Lfj (V) / 

where Z£(V), L^{T) and I^-(V) represent the vis- 
cous terms of the momentum equations, the conduc- 
tion terms of the total energy equation and the viscous 
dissipation term of the total energy equation, respec- 
tively. V = {u,v) is the velocity vector. and 

Lfj are approximated in conservative fashion using cen- 
tral differencing on a standard nine point stencil. In the 
interest of brevity, these are not detailed here. Note, 
for constant viscosity and conductivity, the operators 
and L are linear. 

To facilitate advancing the solution in time, Eq. 
(3) is expressed as a two step scheme as follows: 

!/*• = U?j ~ —Ltj(U n ) + ^ r L v ij (U n ) (8a) 

CFij 

£tf +1 - ^-Z£(I/ n+1 ) = Utf (86) 

The first step, namely the evaluation of E/£, is a purely 
explicit computation. The second step, involving the 
implicit part of the viscous operator, requires resolution 
of a system of equations. This step can be broken up 
into smaller equation sets as follows: Since the viscous 
terms do not contribute to the continuity equation, Eq. 
(8b) can immediately be cast as 


77 n + 1 / 2 

U i±l/2J } S 


rrn. *idU,&tdU 

v 2 2 dt 



(9a) 


Using Eq. (1), we can write 


j-rn+l/2 

U t±l/2,j,S 


' } 2 


( 6 ) 


_&t_(dFy_ dF” dF[ dF*\ 

2 Jij V + dr] dq ) 

In the above, the subscript S stands for “L” when state 
at face (i - 1/2,;) is evaluated and “R” when state at 
face (* + 1/2,;) is evaluated. The complete formula for 
• c is then obtained by appropriately discretizing 
the derivative terms in Eq. (6). For the viscous terms, 
LVj(U n ) (described below) is used unchanged, whereas 
the convection terms are treated exactly as shown in 
Ref. 14. 


V n +1 ^ rm/ V n +i\ _ — V 

v<i 2 ij{ )_ pr 


(96) 


(pe)” +1 - (T?,(T n+1 ) + 4(V n+1 )) = (pc)y 

(9c) 

If viscosity is constant, or if it can be lagged in time 
behind the conserved variables, Eq. (9b) effectively de- 
couples from Eq. (9c). It can thus be inverted first. 
Once Eq. (9b) has been inverted and 1 is known, 
Eq. (2) can be used to recast Eq. (9c) as 




Lij(T n+1 ) = T£ (10) 
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where 


T-i = (7 - 1) {^pt - \v?+ l • vr + r<;(V n+1 ) j 

( 11 ) 

Both Eq. (9b) and Eq. (10) are linear systems 
of equations. They can be inverted using a variety of 
methods. Here, a Gauss-Seidel relaxation scheme is 
used with red-black ordering of the unknowns. Multi- 
grid scheme can be used very effectively to accelerate 
convergence of this relaxation scheme. 

AMR Algorithm for 
Structured Curvilinear Grid Systems 

One of the cornerstones of the present methodol- 
ogy is the use of solution-adaptive mesh refinement in 
structured, body-fitted grid systems. In this approach, 
a relatively uniform body-fitted grid system is used to 
provide the initial tessellation of the flow domain. Since 
solution-adaptive refinement is to be used, there is little 
need during the grid generation to cluster grid points 
near walls or in other regions where high resolution of 
the flow field is anticipated to be needed. Consequently, 
the structured grid system can be optimized with re- 
spect to smoothness and orthogonality alone. In gen- 
eral, high quality grids can be obtained this way since 
grid spacing is not a primary concern during the grid 
generation. 

The AMR Algorithm 

The grid refinement algorithm that best comple- 
ments the use of structured grid systems is the AMR 
algorithm of Berger and Colella. 13 In the AMR algo- 
rithm, the computed solution exists on a sequence of 
mesh levels with finer and finer grids (see Fig. 1). The 
coarsest mesh level covers the entire physical domain 
of interest. Each finer mesh level is created by refining 
cells on the next coarser level in regions where greater 
resolution of the flow field is judged to be needed. A 
fine mesh levels typically covers only a small part of 
the domain and is contained in its entirety within the 
next coarser level. Furthermore, the boundary of the 
fine mesh level must lie a certain distance (measured in 
number of cells) from the boundary of the next coarser 
level, except where the boundaries of both levels coin- 
cide with the boundary of the physical domain. This is 
called “proper nesting” of the mesh levels. Proper nest- 
ing is required so that sufficiently accurate boundary 
conditions can be supplied to the fine grid at coarse-fine 
grid interfaces (see e.g., Ref 13,17). Each mesh level is 
a union of topolgically rectangular blocks. Therein lies 


the congruence between the use of structured grid sys- 
tems and the AMR algorithm. This allows the block 
to be the unit of operation in a computer code. Some 
flexible data structures are needed in the codes to allow 
blocks to be created and deleted. These data structures 
require only small overhead in the codes since the num- 
ber of blocks is relatively small. Within each block, 
however, regular arrays can be used very effectively, 
contributing to the overall efficiency of the method. 

In the AMR algorithm, the generation and main- 
tenance of the hierarchy of mesh level takes place es- 
sentially as follows: Assume there exist a hierarchy of 
mesh levels and solutions on those levels (in the begin- 
ning, the hierarchy is simply the starting coarsest grid 
and the solution is the initial condition). At regular 
intervals, error in each cell on a given level is estimated 
using some suitable measure. Cells where the error ex- 
ceeds a specified limit are tagged for refinement. Errors 
on all finer levels are also checked and cells tagged for 
refinement. Next, tags on the finer levels are propa- 
gated “down” to the subsequently coarser and coarser 
levels, until the level where the error estimation started 
is reached. During this process, additional cells are 
tagged on the coarser levels to ensure that any finer lev- 
els will be properly nested. Once the tagging of cells is 
completed, a special algorithm (see Ref. 13) fits a lim- 
ited number of topolgically rectangular blocks around 
the tagged cells on each level. The mesh within these 
blocks is then refined by an integer ratio to create new 
fine mesh levels. The solution in the old fine mesh levels 
is then copied onto the new mesh levels where the two 
overlap. In regions where new fine grid cells have been 
created, the fine solution is obtained by interpolation 
on the next coarser grid. Note, by this procedure, fine 
grid cells can be removed as well as created. Note also 
that the mesh level where the error estimation started 
does not change. In this work, error is estimated us- 
ing a procedure based on Richardson extrapolation (see 
Ref. 13). 

Refinement of Curvilinear Grids 

In the current methodology, refined curvilinear 
grid systems must be generated from an existing coarser 
curvilinear grid. This grid refinement must be done 
carefully to ensure sufficiently smooth grids on all levels 
of refinement. It is not sufficient to simply connect grid 
points on the coarse grid by straight lines and divide 
the lines into equal segments. This would lead to non- 
smooth grids on the finer levels and to loss in accuracy. 
Therefore a more elaborate grid refinement scheme is 
needed to ensure sufficient smoothness of the fine grids. 
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In this study, the grid refinement is done by com- 
bining parametric cubic spline interpolation and sim- 
plified Hermite interpolation. The cubic splines, here 
natural cubic splines, are used to “reconstruct” the grid 
lines from the discrete grid points. The interpolation 
is done grid line by grid line and produces cubic poly- 
nomials that bridge between any two neighboring grid 
points on a grid line. The parameters used in con- 
structing the splines are simply the coordinates of the 
computational domain. This choice automatically cap- 
tures any nonuniformity of the grid spacing (e.g., grid 
stretching) in the original grid. Thus, when new grid 
points are needed on the grid lines of the original grid, 
uniform spacing in the parameter produces a smooth 
grid. When new grid points are needed that lie in the 
interior of a coarse grid cell, the simplified Hermite in- 
terpolation is used to bridge between the four cubic 
polynomials that define the edges of the coarse grid 
cell. Full Hermite interpolation allows the enforcement 
of both the shape of the edges and also direction of grid 
lines transverse to the edges (transverse derivatives). In 
this usage, however, it is not necessary to enforce trans- 
verse derivatives. Hence the simplified version is used. 
As in the construction of the cubic splines, the param- 
eters used in the Hermite interpolation are the compu- 
tational coordinates of the coarse grid system. Uniform 
spacing in these parameters produces a smooth refined 
grid in the interior of the cell. Note, since the poly- 
nomials describing the shape of the edges of the cell 
were constructed using cubic splines, the overall refined 
grid system, obtained by refining the coarse grid cell by 
cell, will be smooth and at least C 1 continuous. Math- 
ematical formulation of the grid refinement scheme is 
as follows: 

Assume we want to refine cell (i,j). The corners 
of the cell are the nodes (i, j), (i + hj)y {*,j + 1)> an d 
(i + 1, j + 1). The physical coordinates of node (i, j) is 
denoted by r^. The tangential derivative in the direc- 
tion of the i-coordinate is denoted by Sjj, whereas the 
tangential derivative in the direction of the j-coordinate 
is denoted by t^. The and t ij are obtained through 
the cubic spline interpolation. Let £ be the parameter 
that varies along the i-grid lines and 77 be the param- 
eter that varies along the j grid lines. Without loss of 
generality, we can assume that (£,77) = (0,0) at node 
(i, j) and (£,77) = (1, 1) at node (i + 1, j + 1). Now, the 
parametric formulas for the edges of cell (i, j) can be 
written as 


r »,j - i/ 2(£) = ryMO +r»+i,iMO 

+ 81/^13(0 + s *+i.i ^4(0 (12a) 


r »,i+i/2(0 = itj+iMO + r «+i,i+iM0 


+ St,j+1^3(0 + St+l,j+lA-4(£) 

(126) 

r i-i/ 2 j{ 7 }) = + rij+iM 7 ?) 


+ Ujhz{r)) + tij+ih 4 {ri) 

(12c) 

r i+i/2,i( t ?) = r i+ lj'M 7 ?) + rt+i.j+iM 7 ?) 


4- + *1+10+1^4(7?) 

(12d) 

where 


h\(u) = 1 — 3 u 2 -1- 2 u z 

(13a) 

h 2 (u) = 3 u 2 — 2 u 3 

(136) 

hs(u) = u — 2u 2 + u 3 

(13c) 

h±(u) = ~u 2 + u 3 

(13d) 


are the Hermite interpolants. Note, 

M o) = i, h 1 (i) = h\(0) = h' 1 (l) = 0; 

Ml) = l, M0) = *a(0) = Aa(l) = 0; 

M o) = i, MO) = Ml) = *4(1) = 0; 

K(i) = i, mo) = Mi) = >4(0) = 0 . 

Now simplified Hermite interpolation can be written as 

r(£,0 = rjj-xMOMO + r «,j+ i/ 2(OMO 

+ r«-i/2, iOlJMO + r t+i/2,i(’?)> l 2(0 ( 14 ) 

- r^MOMO - r i+i,jM OM 7 ?) 

- r^j+iMOM 7 ?) - r.+W+iMOM’?) 

The function r(£,7j) effectively reconstructs a transfor- 
mation from the computational domain to the physical 
domain. Using this transformation, the grid points for 
the fine grid are placed at uniform intervals in the pa- 
rameters £ and 77. 

The method described above for refining the grid 
system produces smooth grids as it was designed to. 
However, a small complication arises if the flow solver 
is designed to work with cells whose edges are taken to 
be straight lines. In this case, the fine grid cells “seen” 
by the flow solver do not in general match with the 
coarse grid cells. This is illustrated in Fig. 2. This 
complicates any inter-grid transfer operators used in 
interpolation of coarse grid data to the fine grid and 
in coarsening of the fine grid data, particularly if con- 
servation mass, momentum and total energy in the op- 
eration is to be ensured. In Ref. 18, efficient transfer 
operators were designed for this purpose. An alterna- 
tive approach used here is to let the flow solver take 
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the shape of the cells into account. For 2-D problems, 
this only requires that the area of the cells be computed 
using the correct shape. For the formulation used here 
to refine the cells (see Eq. 12-14), the area, <r, of a cell 
is efficiently computed as follows: 


a - JJ dxdy = 1 JJ div{v)dxdy 


r • uds 


(15) 


where r = (x,t/), C denotes the cell, SC denotes the 
edges of the cell, n is the outward pointing normal to 
the cell. The summation is over the faces of the cell. 
The direction of integration around the cell is taken to 
be counter clockwise. This makes the normal n point to 
the right relative to the direction of integration. After 
a little algebra it can be shown that when the cells are 
defined by the cubic polynomials, the integral over a 
cell face is 


- / r • nds = ~{x a yb ~ XbVa) 

1 J 6C\ 1 

+/ (r OI (dr/ds) a ,T b , (dr/ds) b ) ( 16 ) 

where 

/ (r„, (dr/ds) a ,r b , ( dr/ds ) b ) = 

60 \\ds) a \ds) b \ds) h \dsJJ 

X-(S)J) ™ 

-35 (<»-*•) ((£),- (s)J) 

In the above, subscripts “a” and u b” indicate the begin- 
ning node and end node on the edge, corresponding to 
the direction in which the parameter “s” in the integral 
varies. Note, the first term in Eq. (15) corresponds to 
integration along a straight line from r a to r^. The sec- 
ond term, therefore, can be thought of as a correction 
to the integral corresponding to the deviation of the 
curve from a straight line curve. Changing the direc- 
tion of integration in Eq. (16) (while keeping a normal 
pointing to the right) simply changes the sign of the 
results, i.e., 

f (r a , ( dr/ds) a , r b , (dr/ds) h ) = (18) 

- / (r», - ( dr/ds) b , r a , - (dr/ ds) a ) 


Applying Eq. (14-16) to cell (t, j), and keeping 
with the notation of Eq. (12-14), we obtain 

*ij = °ij 

“b / ( r »+i,i>^*+iii> r «+iii+ii^*+iii+i) (1®) 

~b / ( r i+l,j+lj “Si+ljj+l? r i,j+l> 

■b / r *,ii “ ^t,y) 

where 

a \j = 2 “ x *,j) (y*+i,i+i Vij) (20) 

is the area of a cell whose faces are straight lines. Equa- 
tion (19) can finally be written as 


where 


(Tij 




+ (/t+l/2,j “ /t-l/2,i) 

- (/ »,jf+l /2 “ fij-l/ 2 ) 


( 21 ) 


fi+ 1 / 2 , j “ / j^i+iii+i) (22a) 

fij+ 1/2 = / (226) 

As Eq. (21) and (22) suggest, a ij for a block of cells 
is most efficiently computed by evaluating /i+ 1 / 2 ,; and 
fi,j+ 1/2 for the edges of the mesh and adding the cor- 
rection to cr*j . 

Note, to ensure that cell edges on the coarsest grid 
level always coincide with grid lines on the finer levels, 
the grids on all finer levels are obtained by refining 
the coarsest level. Consequently, the spline coefficients 
used to reconstruct the grid lines need only be known 
for the nodes on the coarsest level. 


Transfer of Data Between Levels 

When new fine grids are created the solution on 
that grid must be initialized by interpolating the data 
on the underlying coarse grid. Interpolation from coarse 
grids is also needed at interfaces between coarse and 
fine grids to provide boundary conditions for the fine 
grid. Also, when fine grids are deleted, the data from 
those grids must be transferred to the underlying coarse 
grid. In all cases the transfer of data must ensure con- 
servation of mass, momentum and total energy in or- 
der to maintain accuracy and, for example, to maintain 
correct speed of moving discontinuities such as shocks. 
With the "current definition of cells, the transfer of data 
from a fine grid to a coarse grid is accomplished by sim- 
ply adding up the mass, momentum and total energy 
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in all the fine grid cells that correspond to a particu- 
lar coarse grid cell. Transfer of data from coarse to fine 
grids is currently done using a conservative linear inter- 
polation as shown in Ref. 18. That particular approach 
works well on relatively uniform grids but may need to 
be improved for grids with rapidly changing cell areas 
(Jacobians). 

Implementation — Object-Oriented 
Mixed Language Programming 

Implementation of a solution adaptive mesh re- 
finement algorithm like the AMR algorithm described 
above, requires the use of programming languages that 
support dynamic memory allocation (and de-allocation) 
and user-defined data structures. The former capabil- 
ity is needed so that mesh levels and blocks can be 
created and deleted extemporaneously as the solution 
developes in time. The latter capability is desirable so 
that effective organization of the data can be done sys- 
tematically and autonomously in the computer codes. 
Programming languages that offer both capabilities in- 
clude C and C++. FORTRAN90 will also offer some 
of those capabilities. 

While dynamic memory management and flexible 
data structures are needed for an effective implemen- 
tation of the methodology, efficient floating point op- 
erations are also needed for fast execution of the code. 
The programming language that currently offers the 
most efficient floating point operations, particularly on 
vector supercomputers, is FORTRAN77 (due to highly 
developed compilers). This language does not, how- 
ever, have the needed capabilities for memory manage- 
ment and data structures. Fortunately, the modularity 
of the AMR algorithm allows one to take advantage of 
the strength of the different programming languages. 
FORTRAN can be used very efficiently to implement all 
operations within a block that are related to advancing 
the solution in time, computing fluxes, applying physi- 
cal boundary conditions, etc. A driver module that al- 
locates memory for blocks, controls the time stepping, 
error estimation and refinement, and calls the FOR- 
TRAN routines can then be implemented in another 
programming language. 

In this work, the AMR driver module was imple- 
mented using the C++ programming language. This 
language is very well suited for this purpose due to its 
support for object oriented programming and well de- 
fined procedure for calling FORTRAN programs. As 
the AMR algorithm suggests, the basic object in the 
implementation is the block. The corresponding class 
in the code is called “LevelBlock.” A class in C++ is a 


user defined type and consists of data and a collection 
of functions (member functions) that operate on the 
data. In this case, the data in “LevelBlock” consist of 
arrays for the conserved variables, primitive variables 
and metrics, and a special data type for information 
about physical boundary conditions for the block. The 
member functions include functions that call the FOR- 
TRAN implemented flow solver. A second object in the 
implementation is a class called “MeshLevel.” The data 
in this class includes the collection of “LevelBlocks” 
that make up a single mesh level, and pointers to the 
“MeshLevel” objects that contain the next coarser and 
next finer mesh levels. In this work, extensive use was 
made of the AMR library developed by Crutchfield and 
Welcome. 19 The AMR library is a collection of dimen- 
sion independent classes specially designed to aid in 
implementation of schemes employing the AMR algo- 
rithm. 

Results 

When this paper is written, the methodology de- 
scribed in previous sections has only been tested on 
a number of test cases involving inviscid flow. Here, 
three such cases are presented. These are a subsonic 
flow over a NACA0012 airfoil, a transonic flow over a 
NACA0012 airfoil, and a supersonic flow over a blunt 
body. In all cases only one level of refinement is used 
with a refinement ratio of four. 

NACA0012 airfoil at M = 0.5 

The first test case is a subsonic flow over a NAC A- 
0012 airfoil at zero angle of attack. The free stream 
Mach number (M) in this case was taken to be 0.5. To 
take advantage of the symmetry of the geometry, the 
computations were restricted to the half-plane above 
the airfoil. The starting (coarse) grid used in the com- 
putations contained 32 by 72 cells. The far field bound- 
aries of the grid system were between 50 and 100 chords 
from the airfoil. Figure 3 shows the center region of the 
grid to a distance of about 4 chord lengths from the air- 
foil. This coarse grid has only 12 cells over the surface 
of the airfoil (on one side), and 48 cells after refinement. 
Figure 4 shows the region around the airfoil that was 
refined and Fig. 5 shows contours of pressure coefficient 
near the air foil. The refined region consisted of three 
blocks. A total of 9.5% of the coarse grid cells were 
refined. The computed solution compares well with ex- 
perimental data published in Ref. 20. According to 
this data the minimum pressure coefficient at the sur- 
face is -0.4687, compared to computed value of -0.471, 
an error of about one half a percent. 
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NACA0012 airfoil at M = 0.8 

The second test case is a transonic flow over a 
NACA0012 airfoil at zero angle of attack. The free 
stream Mach number was taken to be 0.8. The start- 
ing grid system was the same as described above and 
shown in Fig. 3. Figure 6 shows the refined grid system, 
and the computed pressure coefficient is shown in Fig’s 
7 and 8. As seen in Fig. 6, the refined grid consisted 
of two blocks. Fewer than 12% of the coarse grid cells 
were refined, corresponding to a savings by a factor of 
more than eight if the entire grid had been refined. The 
present scheme resolves the shock on the airfoil very 
crisply, within only two cells. As seen in Fig’s 7 and 8, 
the location of the shock is at about 48% chord. Ac- 
chording to experimental data reported in Ref. 20, the 
correct shock location is 40% to 44% chord. The rea- 
son for the discrepancy is lack of resolution — experience 
with similar geometries indicate that with one extra 
level of refinement located right over the shock, the 
correct shock location and strength will be predicted 
with the current scheme. 


Supersonic flow over a blunt body at M — 5 

The last test case is a supersonic flow over a 
26.56° wedge with a round leading edge. The lead- 
ing edge is a cylindrical surface whose radius is 0.125 
in the current scale. The free stream flow is at M = 5 
and a detached bow shock is formed. 

The grid system at the end of the computations 
is shown in Fig. 9. Figures 10 and 11 show contours 
of density and entropy, respectively, around the wedge, 
whereas, Fig. 12 shows plot of density on the stagna- 
tion streamline versus axial location. As Fig. 9 shows, 
the refined grid consisted of several topologically rect- 
angular boxes. The refinement has mainly taken place 
around the shock and the stagnation region. Approxi- 
mately 47.6% of the coarse grid cells were refined. As 
Fig. 11 to 12 show, the shock is very well captured. 
Also, the entropy generated at the shock is properly 
convected along the stream lines. In Fig. 12 it can be 
seen that no overshoots or undershoots in density are 
formed at the shock, and the density increases mono- 
tonically behind the shock. As Fig. 12 shows, the shock 
stands at axial distance between 0.060 and 0.063 from 
the leading edge. This corresponds to ratio of stand-off 
distance over two times the radius of curvature at the 
leading edge of 0.241 to 0.252. In comparison, Ref. 21 
shows an experimental value of about 0.24 for flow at 
M — 5 over a cylinder. 


Conclusions 

In this paper, a methodology is proposed for sim- 
ulating unsteady viscous and inviscid flows. The cor- 
nerstones of the methodology are the use of structured 
(multiblock) grid systems, solution adaptive mesh re- 
finement based on the AMR algorithm of Berger and 
Colella 13 , and a hybrid explicit-implicit discretization 
of the Navier-Stokes equations. The methodology is im- 
plemented in a computer code which is written using 
mixed language programming — C++ for a driver mod- 
ule and FORTRAN for a flow solver module. At this 
point only limited tests of the methodology have been 
performed. Overall, the results of those tests were very 
good. However, it was found that the simple interpola- 
tion scheme used to transfer data from coarse grids to 
fine grids did not work sufficiently well on meshes with 
rapidly changing grid spacing in two directions at the 
same time. Further development in this area is needed. 
The development and testing of the methodology is on- 
going. 
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Figure 1. A grid system on three levels — a coarse 
grid covering the entire physical domain and two prop- 
erly nested fine grid levels, each consisting of several 
topologically rectangular blocks. 



Figure 3. Initial (coarse) grid system for NACA0012 
airfoil — computations were done on only the upper half 
of the symmetric grid. 



Figure 2. A curvilinear coarse grid system and a re- 
fined grid — cells are defined by straight lines between 
nodes, resulting in fine grid cells that do not match 
with the “parent” cells on the coarse grid. 
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Figure 4. Refined grid around a NACA0012 airfoil — 
flow at M — 0.5 and zero angle of attack. 
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Figure 5. Flow over a NACA0012 airfoil atM = 0.5 
and zero angle of attack — contours of pressure coeffi- 
cient. 



Figure 0. Refined grid around a NACA0012 airfoil — 
flow at M = 0.8 and zero angle of attack. 
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